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Abstract 



We present a perturbation theory for non-Markovian quantum state diffusion 
(QSD), the theory of diffusive quantum trajectories for open systems in a bosonic 
^ I environment [Physical Review A 58, 1699, (1998)]. We establish a systematic ex- 

^^' pansion in the ratio between the environmental correlation time and the typical 

system time scale. The leading order recovers the Markov theory, so here we con- 
k> , centrate on the next-order correction corresponding to first-order non-Markovian 

^ ' master equations. These perturbative equations greatly simplify the general non- 

Markovian QSD approach, and allow for efficient numerical simulations beyond the 
Markov approximation. Furthermore, we show that each perturbative scheme for 
QSD naturally gives rise to a perturbative scheme for the master equation which we 
study in some detail. Analytical and numerical examples are presented, including 
the quantum Brownian motion model. 
PACS Numbers: 03.65.-w, 03.65.Bz, 05.40. -hj, 42.50.Lc 



1 Introduction 

Recently, a non-Markovian quantum trajectory theory - named non-Markovian quantum 
state diffusion(QSD) - that describes the dynamics of a quantum 'system' coupled to an 
'environment' of harmonic oscillators has been presented [1]. Many outstanding new ex- 
perimental advances can be properly studied only if non-Markovian effects are taken into 
account. These include experiments with high-Q microwave cavities, quantum optics in 
materials with a photonic bandgap, or output coupling from a Bose-Einstein condensate 
to create an atom laser [2, 3, 4, 5, 6], to name a few. Also, the important phenomenon 
of decoherence which takes place on time scales that can be of the same order as the 
correlation time of the environment require theories beyond the standard Markov approx- 
imation. Further motivation are more fundamental questions about the proper description 
of individual open systems in quantum mechanics. Indeed the infamous problem of the 
"Heisenberg cut" (understood here as the cut between the system and the environment) 
is intimately related to the non-Markovian evolution of the system when the environment 
is ignored. 

In the Markov regime, quantum trajectory approaches using stochastic Schrodinger 
equations for pure states of the system play an important role in quantum optics, par- 
ticularly for numerical simulations [7, 8, 9, 10, 11, 12, 13]. These Markov models also 
have many appealing features from a theoretical and conceptual point of view [8, 9, 10, 
14, 15, 16]. It is therefore desirable to generalize the powerful quantum trajectory ap- 
proach from the Markov regime to the more general case of non-Markovian evolution. 
Several attempts have been made recently from different perspectives. The linear non- 
Markovian unravelling underlying our theory was developed in [18, 19, 20] - see also [17] 
for a related attempt. Alternatively, a non-Markovian theory based on pseudomodes and 
a non-Markovian jump approach have been developed recently [21, 22, 23, 24, 25]. 

In ref. [1], the ultimate nonlinear non-Markovian QSD equation relevant to this paper 
has been derived directly from a microscopic model. In this framework, the reduced 
density matrix of the subsystem obtained by tracing over the environmental degrees of 
freeedom is unraveled into an ensemble of continuous trajectories which correspond to 
the various realizations of the driving complex Gaussian process. Non-Markovian QSD 
has been applied to interesting and physically relevant models where both computational 
power and many new features have been demonstrated [1]. However, many issues are still 
to be addressed. In particular, further applications of non-Markovian QSD to a variety 
of realistic problems are desirable. In addition, the theoretical implications of this new 
approach remain to be explored. Clearly, a full exploration of non-Markovian QSD is 
an extensive project. The purpose of the present paper is a step towards this extensive 
project. 

So far, although being a general theory, it is difficult to implement non-Markovian QSD 
directly on a computer in all generality. It has been applied to a variety of model problems 
in [1]. However, in this paper we show how the non-Markovian QSD approach allows to 
find a systematic expansion of the reduced system dynamics in powers of the ratio between 
the environmental correlation time and typical system time scales. Thus, in order for non- 
Markovian QSD to have more applications, we establish a useful and practically relevant 
perturbative approach that is directly amenable to computer simulations. Our first and 



most important motivation for this paper is to present a systematic perturbative approach 
for non-Markovian QSD around the Markov hmit. This perturbative "post-Markovian" 
QSD scheme is a time-dependent approach which preserves the non-Markovian nature of 
the problem in each order of approximation. 

The second purpose of the present paper is to estabhsh the relationship between the 
non-Markovian QSD equation and non-Markovian master equations. Such a relation is 
well-known in the Markov regime, where one generally starts from the standard Lind- 
blad Markov master equation to read off the Markov QSD equation [9, 10]. For non- 
Markovian dynamics, closed master equations are rare and thus, the direct link between 
non-Markovian QSD and the corresponding master equation is of great interest. In order 
to obtain the master equation from its non-Markovian QSD counterpart, one has to take 
the ensemble mean over the stochastic process driving the QSD trajectories analytically 
(see also ref. [1]). In general, this is only possible for simple models. In our pertubative 
scheme, however, we are able to derive master equations directly from non-Markovian 
QSD which turns out to be useful from both theoretical and practical points of view. It 
is important to note that a master equation derived in this way will necessarily preserve 
the positivity because such density matrices can be decomposed into pure states at all 
times. We thus also address the problem how to ensure positivity for non-Markovian 
master equations - a difficult subject in itself. 

The third motivation of the paper, therefore, is to present a perturbation approach to 
non-Markovian master equations. Using an example we show that the resultant approxi- 
mate master equation preserves positivity. We also analyse the approximations leading to 
the Caldeira-Leggett model [26] - which is known to violate positivity for certain initial 
conditions on short time scales. 

The organization of this paper is as follows. In Section II we briefly review the basic 
concepts and equations of non-Markovian QSD. The aim of this section is two-fold: first, 
to act as a brief introduction for readers not familiar with the subject, second, to serve as 
the natural starting point of our new development. In Section III we develop the formal 
time-dependent perturbation theory for non-Markovian QSD. In Section IV, we present a 
systematic method of deriving the master equation from non-Markovian QSD. We show 
that a perturbative QSD scheme naturally leads to a perturbative scheme for the master 
equation. We will apply the approximation schemes developed in this paper to some 
examples in Section V. In Section VI we take quantum Brownian motion as a typical 
model to illustrate the perturbative schemes based on QSD for the master equation. 
We conclude the paper in Section VII, while some useful material can be found in the 
Appendices. 

2 Non-Markovian quantum state diffusion 

Both Markov and non-Markovian QSD are based on two related stochastic dynamical 
equations, a linear one and a nonlinear one. The linear one is mathematically simpler. 
However, it does not preserve the norm of the state vector, which in general tends to zero. 
Hence, only the nonlinear equation, which preserves the norm, can be interpreted as a 
distribution of time-dependent pure states with given probabilities (ie as an unraveling): 



the density matrix is then given by the ensemble mean of the pure states, at all times. 
Moreover only this nonlinear equation is suitable for numerical simulation [1, 20]. Never- 
theless, we start our presentation with the simpler linear equation, leaving the nonlinear 
one for the following subsection. 

2.1 Linear non-Mar kovian quantum state diffusion 

Our quantum trajectory theory is based on a standard model of open system dynamics: 
a quantum system interacting with a bosonic environment with total Hamiltonian 

i^tot = H + Y,g\{La\ + L^aa) + JZ^AflW, (2.1) 

A A 

where H is the Hamiltonian of the system of interest and L, a system operator coupling to 
environment, is called here Lindbald operator (as it plays the role of a 'Lindblad operator' 
in the Markov limit). The linear non-Markovian QSD equation [18, 19, 20] unravelling 
the reduced dynamics of model (2.1) takes the form^ 

^^^ = -tHiJt + LztiJt - Lt /* a(t, s)^rfs, (2.2) 

at JO ozs 

where Zt is a colored complex Gaussian process with zero mean and correlations 

M[z;z,] = a{t, s), M[ztZs] = 0. (2.3) 

The bath correlation function a{t,s) in (2.2) has to be a positive Hermitian kernel: 
a{t,s) = a{s,t)*. This non-Markovian unraveling is ensured to be consistent with the 
reduced density operator approach since the ensemble mean of the solutions of Eq. (2.2) 
over the noise Zt will reproduce the density matrix of the system, 

pt = Trenv (e-^^-*|V^o)(V^o| ® pre^""'"*) = M [|^i(^))(Vt(^)|] . (2.4) 

Here M[- ■ ■] denotes the ensemble average over the classical driving noise Zf. 

From Eq. (2.2), we see clearly that the evolution of the state ipt at t depends on the 
whole history of the noise z. The equation (2.2) can be written in the more appealing 
form 

— ^t = -iH^pt + Lzt^pt - L'^ I Oi{t, s)0{t, s, z)ds^t (2.5) 

at Jo 

by defining an operator ^ 0{t, s, z) in (2.5) such that 

d{t,s,z)ij,^^^. (2.6) 

It turns out that 0{s,s,z) = L. The t-dependence of the operator 0{t,s,z) can be 
determined by the consistency condition 

d 6 _ 5 d 
dt 6zs Szs dt 



^For simplicity, wc choose ft = 1 throughout the paper. 

^ The notation 0(i, s, z), rather than 0{t, s, zt) reflects that the operator 0(i, s, z) contains the noise 
z in a nonlocal way, that is, it might be dependent on the whole histories of the noise {zg : Q < s <t\. 



together with the hnear non-Markovian QSD equation (2.5). 

The appeal of equation (2.5) (or (2.2)) is its hnearity which is often useful in the 
mathematical analysis of our QSD approach (see Section IV). Its use as simulation tool 
is severely undermined since, for an infinite heat bath, the norm \\ipt\\ of the solutions 
of (2.2) goes to zero with probability one and to infinity with probability zero. For this 
reason, an unraveling in terms of normalized states is crucial for non-Markovian QSD to 
be truly useful for numerical simulations. 

2.2 Nonlinear non-Markovian quantum state diffusion 

The non-Markovian QSD unraveling based on normalized states 

M^) = ^ (2.8) 

has been derived recently [1] from the linear non-Markovian QSD equation (2.5) by making 
use of a Girsanov transformation of the noise. We get 

-Jjt= - iH^t + {L-{L)t)^tZt 
at 

'a{t,s)\{L^-{L^)t)d{t,s,z) 

L 

- {{L^-{L^)t)d{t,s,z))t\ds4^t (2.9) 



where Zt is the shifted noise. 



zt = zt+ f'a{t,sy{L^)Js, (2.10) 



and {L)t = {ipt\L\ipt) denotes the quantum average. Again here, the ensemble average of 
the solution to Eq. (2.9) recovers the density matrix of the system, 

Pt = M[\i,t{z)){Mz)\\. (2.11) 

This nonlinear non-Markovian QSD equation can be rewritten in a more compact form: 

^^i = -tH^t + ML)^tZt - ML^)0{t, z)^t + {ML^)0{t, ~z))tk (2.12) 

at 

where A((A) = A — (A)^ for any operator A and 

0{t,z)= / ci{t,s)0{t,s,z)ds. (2.13) 



Equation (2.9) (or (2.12)) is the fundamental equation of non-Markovian QSD ^, and 
is our starting point for the perturbative approach in Section III. Numerical simulations 
of non-Markovian open system dynamics using this equation can be found in ref. [1]. 

■^From now on, unless otherwise emphasized, non-Markovian QSD refers to equation (2.9) or (2.12), 
not to the hnear equation (2.2). 



Note that the non-Markovian QSD equation (2.9) reduced to the standard Markov QSD 
equation [10] for a(t, s) -^ 6(t — s). 

To conclude this section, we would like to make two remarks about the non-Markovian 
QSD approach. First, the derivation of both the linear and non-linear non-Markovian QSD 
equations are based on the assumptions that the environment is bosonic, and that initially 
the state of the total system+environment is factorable Po = Po ® /^o"*') where the initial 
state of the system p^ = \i/jo) {ipol is independent of the noise Zf. In fact, if L 7^ L^, (2.9) 
is valid for a zero temperature environment only and the equation for finite temperature 
gets two additional terms, see [1]. Secondly, since Eq. (2.2) and hence (2.9) are derived 
directly from the microscopic model, these QSD equations can be read off automatically 
from the total Hamiltonian (2.1). This suggests that the non-Markovian QSD approach 
also represents a brand new way to derive the quantum master equation of open quantum 
systems (see Section IV). 

3 Time-dependent perturbation theory 

The non-Markovian QSD approach offers a very promising method to handle quantum 
systems whenever non-Markovian effects are relevant. However, many interesting prob- 
lems which arise in open quantum systems are such that the operator 0{t, s, z) appearing 
in Eq (2.9) cannot be determined exactly. Moreover, the nonlocal noise contained in the 
fundamental equation (2.12) might cause difficulties in numerical simulations. In this 
section, we aim for a formal perturbation scheme for the non-Markovian QSD equation 
(2.9). Applications of the general perturbative method developed here will be presented 
in Section V. 

3.1 First order approximation of the operator 0{t^s^z) 

We need to know the operator 0(t, s, z) from (2.6) in order to solve the non-Markovian 
QSD equation (2.9) on a computer. Notice that O enters (2.9) under the memory integral 
/q Q;(t, s)0(t, s, z) ds only. Therefore, if the correlation time of the environment is not too 
long, only s- values in the vicinity of the upper integration range are relevant, and thus 
we consider the expansion of the operator 0(t, s, z) in Eq.(2.9) in powers of (t — s). 



A/ X A/ X dO(t,s,z) 

0{t,s,z) = 0{s,s,z)+ • ^ ' ' ^ 



+ 



dt 

l d'^d{t,s,z) 
2 dT^ 



(t-s) 



t=s 
2 



it-sy + ... (3.1) 



t=s 



Substituting (3.1) into (2.9) or (2.12), we get a hierarchy of approximate QSD equa- 
tions by truncating the above expansion. The validity of the corresponding approximation 
depends on the environment correlation time r determined by the correlation function 
a{t, s). In fact, it turns out that O changes on system time scales as a function of (t — s) 
and thus, the expansion (3.1) corresponds to a systematic expansion of the non-Markovian 
QSD equation in powers of the number ujt where a; is a typical 'system' frequency and r 



the environmental correlation time. For example, the zeroth-order term leads to the stan- 
dard Markov QSD, when r -^ 0. The first-order term is the most important correction 
to the Markov dynamics. Therefore, in what follows we will work out the approximation 
up to the first order in some detail. 

By using the consistency condition (2.7), one can work out the following expressions 
for the operator 0{t,s,z) at time point t = s, without knowing its explicit form (for 
details, see Appendix A): 

d{s,s,z) = L (3.2) 

dd{t,s,z) 



dt 



t=s 



-i[H,L]- r a{s,u)du[L\L]L (3.3) 

Jo 



where H is the Hamiltonian and L is the Lindblad operator, as specified in the previous 
section. Now, we are in a position to write out the non-Markovian QSD equation up to 
the first order. Indeed, the first two terms in the expansion (3.1), substituted repectively 
by (3.2) and (3.3), yield 0{t) in (2.13) in the following form: 

0{t)=go{t)L-g^mH,L]-g2{t)[L\L]L (3.4) 

where 

go{t) = / a{t,s)ds (3.5) 



'0 

^i(t) = / a{t, s){t - s)ds (3.6) 

Jo 

92{t) = a{t, s)a{s,u){t ~ s)duds (3.7) 

Jo Jo 

Note that go is of the order one, yet gi and g2 are of the order of the environmental 
correlation time r. Substituting (3.4) into (2.12), the first-order non-Markovian nonlinear 
QSD equation is obtained: 

—Jt= - iHi,t + ^t{L)i,tZt 
dt 

- go{t){At{L^)L-{AtiL^)L)t)^t 

+ igrit) [ML^)[H,L] - {/\t{L^)[H,L\)t)^t 

+ (72(t)(A,(Lt)[Lt,L]L-(A,(Lt)[Lt,L]L),)v^, (3.8) 

where Zt is the shifted noise, At{L) = L — {L)t, and {L)t = {ipt\L\ipt) is the quantum 
expectation value. 

The Hamiltonian H defines a typical system frequency u, the combination L'^L defines 
a typical system relaxation rate F. We thus see that the zeroth order term in (3.4) gives 
rise to a term of the order F (second line in eq.(3.8), whereas the two first order terms in 
(3.4) lead to corrections which are smaller by a factor ur or Fr, respectively (third and 
fourth line in eq.(3.8), r is again the environment correlation time). Therefore, we expect 
(3.8) to be valid for non-Markovian situations where the environmental correlation time 



may be finite but no larger than typical system time scales. The Markov case emerges 
for r ^ 0, where the first order correction becomes negligible and only the zeroth order 
term remains. Then (3.8) reduces to the standard Markov QSD equation for t > 0. 

We also see that non-Markovian properties are encoded in the time dependent coef- 
ficients gi{t) which change on the very fast environmental correlation time scale r. The 
absence of the noise z in the first order expansion (3.4) is remarkable (higher order ex- 
pansions contain the noise). Note that the approximate non-Markovian QSD equation 
(3.8) still preserves the norm of the wavefunction. Eq. (3.8) is the main result of this 
section. As stated before, applicability of non-Markovian QSD lies in the determination 
of the operator 0{t,s,z). As we have already pointed out, the difficulties in handling 
non-Markovian unravelings is often the nonlocal noise z appearing either in the func- 
tional derivative (see (2.2)) or in the integrand operator 0{t,s,z) (see (2.5) and (2.9)). 
We see that the above approximate QSD equation greatly simplifies non-Markovian QSD 
equation (2.12). 

In addition, Eq.(3.8) is explicitly written in terms of the Hamiltonian of the system 
H, the Lindblad operator L and their various commutators. All of these can be obtained 
automatically once the physical model is specified. The only work left is to calculate the 
coefficients gi{t){i = 0, 1, 2) from the environment correlation function a{t, s). 

After working out the formal perturbative QSD equation, it is useful to see the concrete 
form of the coefficients gi{t). For simplicity, we assume here that the system is driven by 
Ornstein-Uhlenbeck noise, characterized by the exponential correlation function 

«(t,s) = 2e-^l*-l (3.9) 

where 7~^ = r defines the finite environmental memory or correlation time. Note that 
this corresponds to a Lorentzian spectrum. In the limit 7 — > oo, the Ornstein-Uhlenbeck 
noise reduces to simple comples white noise: 

a{t,s) = 5{t-s) (3.10) 

In the case of the Ornstein-Uhlenbeck process, the coefficients gi{t) can be easily obtained 
fromEq. (3.5)-(3.7): 

9o{t) = ^(l-e-^O (3.11) 

9i{t) = ^ (1 - e-^* - 7te-^*) (3.12) 

g,it) = ^ (^1 - e-^* - 7te-^* - ^ft V^*) (3.13) 

In the long-time limit t ^ r, we see that the coefficients of the non-Markovian QSD 
(3.8) become constant: go = 1/2, gi = 1/27,(72 = 1/47, which also confirms that go is of 
the order one whereas gi and g2 are of the order of the environmental correlation time 
r = 7~"^. 

In the Markov limit 7 -^ 00, go{t) -^ 1/2 and gi{t),g2{t) ^ for t > and the non- 
Markovian QSD equation (3.8) reduces to the standard Markov QSD equation [10] (Note 



here we write it in the Stratonovich form [27, 28]): 

^^, = -iHi,^ + A,(L)V^, o {zt + (Lt),) - iA,(LtL)Vii (3.14) 

with Zt the standard complex white noise, as expected. 

Our formal perturbation approach can be carried out to any desired order of approx- 
imation (For the details of the second order expansion and the coefficients, see Appendix 
B). It is important to note, however, that the higher order derivatives of 0(t, s, z) at t = s 
may contain the noise z. 

Since the linear non-Markovian QSD equation (2.2) is often simpler to derive the 
corresponding master equation (see Section IV), we also give its first-order approximation: 

i^t = -iHiJt + L^tzt - goit)L^Li^t + igi{t)L^[H, L]^t + 92{t)L'^[L\ L]LiPt (3.15) 

where the coefficients gQ{t) , gi(t) , g2(t) are given by (3.5), (3.6) and (3.7). 

3.2 Functional expansion of 0{t,s,z) 

In this subsection, we consider another kind of perturbative expansion, the functional 
expansion of the operator 0{t, s, z) in terms of noise z^: 

dit,s,z) = do(t,s) 

ft . 
+ / Oi(t,s,v)zydv 
Jo 

+ d2{t,S,Vi,V2)Z^^Z^^dVidV2 



JQ 

+ .... 

+ / ... / dn{t,S,Vi,...,Vn)z^^...Z^^dVi...dVn 

Jo Jo 
+ ... (3.16) 

where the operators 0„(t, s,Vi, ...,f„) are independent of the noise z and are symmetric 
in their n last variables (e.g. 02(t,s,Vi,V2) = 02(t,s,V2,Vi)). The initial condition is 
0(t, t, z) = L. The expansion (3.16) takes into account the generally nonlocal dependence 
of the operator 0{t, s, z) on the noise z. 

From the consistency condition (2.7) and the QSD equation (2.2), we get a hierarchy 
of equations for the operators 0„(t, s,fi, ...,f„) (see Appenix A). Of particular interest 
is the zeroth order term Oo(t, s), which satisfies the following equation(ignoring the first 
order term): 

^Oo(t,s) = [-iH,do{t,s)] 

- [L^Ooit),do{t,s)] (3.17) 

For the approximation 0{t,s,z) ^ Oo(t, s), the approximate non-Markovain QSD equa- 
tion then takes form: 

^Jt = -iHi,t + ^t{LWzt - ML^)Oomt + (A,(Lt)Oo(t))i7/i, (3.18) 

dt 



The justification of this approximation is that whenever the open quantum system deviates 
shghtly from the Markov dynamics, then the first term Oo{t,s) of the expansion (3.16) 
plays the dominant role. This can be easily seen from the fact that all 0„(t, fi, ...,f„) = 
/o a(t, s)0„(t, s,fi, ...,f„)(is,n > 1 go to zero in the Markov limit: a{t,s) -^ 5{t — s), 
except the first term Oo{t) = /q a{t, s)Oo{t, s)ds which goes to |L. Physically this can be 
understood as follows. In the Markov case, the quanta coupled from the system to the 
environment never come back to the system. Whereas, in the non-Markovian case, the 
emitted quanta will re-couple from the environment to the system. 

Similarly, one can build the higher order approximations which usually contain the 
noise z. One obtains then a series of approximate QSD equations. The master equation 
corresponding to the zeroth order approximation (3.18) is derived in the next section. 

4 Non-Markovian QSD versus non-Markovian mas- 
ter equation 

In this section, we discuss how to derive the master equation from the non-Markovian 
QSD equation. Our motivations are as follows. First, the master equation approach has 
a long tradition and is fundamental in open quantum system dynamics, and the reduced 
density operator contains all mean values of the 'system' that can be directly observed 
and measured. Second, although it is clear in principle that each perturbative scheme for 
non-Markovian QSD gives rise to a perturbative scheme for the non-Markovian master 
equation, it is very difficult in practice to carry out this program without a systematic 
way to derive the non-Markovian master equation from its QSD counterpart. The aim 
of this section is to show how to derive the quantum master equation directly from the 
non-Markovian QSD. Based on this result, we establish explicitly the relation between 
the perturbative QSD equations and perturbative master equations. 

4.1 General master equation 

The starting point of the derivation of the general master equation is the unnormalized 
projection operator Pf, 

p, = \Mz))m^)\ (4.1) 

Recall that the reduced density operator can be reproduced by taking the statistical means 
over the noise: pt = M[Pt] = M [\ilJt{z)){ilJt{z)\]. Accordingly, the temporal evolution 
equation for Pf can then be obtained from (2.5): 

j^Pt= - i[H,Pt]+LPtZt + PtLhl 

- L^ f a{t,s)d{t,s,z)dsPt-Pt f a{t,s)*d{t,s,zydsL (4.2) 

Jo Jo 

The above equation is, of course, a stochastic differential equation with time- dependent 
coefficients. Accordingly, the master equation corresponding to Eq. (2.2) may be obtained 
by taking statistical mean values of Eq. (4.2). 



To this end, we note that for any complex Gaussian noise Zt, the following relations 
hold (see Appendix C): 



dsM[ztzl]M 



SPt 



Sz* 

SZs 



M[PtZt] = 

M\Ptzl] = I dsM[zrZs]M 
Jo 

From (4.3) and (4.4), the following identities are obtained 

M[LPtZt] = L I a{t,syM\Ptd{t,s,zy 



(4.3) 
(4.4) 



ds 



M[PtLhn = 

Here we used the following relations: 



a{t, s)M \d{t, s, z)Pt\ dsL^ 



M 
M 



Sz., 



_5_ 
Jz* 



P 



M[—\^t){^t\] = M[d{t,s,z)P,] 



(4.5) 
(4.6) 

(4.7) 
(4.8) 



and we take advantage of the definition of the 0-operator (2.6). The validity of the above 
two identities (4.7) and (4.8) is ensured by the fact that the solution ipt of Eq. (2.2) 
is the analytic function of z and is thus independent of z*. Accordingly, 6\ilJt)/Szl = 
0, 6{^lJt\/Szt = 0. 

Hence, using (4.5) and (4.6), the exact non-Markovian master equation corresponding 
to non-Markovian QSD (2.2) can be obtained: 



1/^ = "^''''^^ + 



L,M 



PtO{t,zy]] -\L\M\0{t,z)Pt 



(4.9) 



Where as before M[- ■ ■] stands for the ensemble average, and 0(t, z) is defined in (2.13). 

Eq. (4.9) is the exact equation on which our perturbation approach is based. As 
an evolution equation, the above master equation does not look very nice since the last 
two terms appearing in the equation have not yet been written in terms of p. It seems 
quite difficult to write this equation into a closed evolution equation in full generality, if 
not impossible. We shall see, however, that in many interesting and physically relevant 
situations, a closed form for this equation can be found (see below). Notably, the use of 
the relations (4.3) and (4.4) can make a tremendous simplification in deriving the master 
equation of open quantum system from its QSD counterpart. In fact, it enables us to 
find out an exact or an approximate non-Markovian master equation by directly using 
the techniques of stochastic process. 

The non-Markovian master equation (4.9), by design, will always preserve the posi- 
tivity, trace and hermiticity. 
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4.2 Approximate master equations 

Since the master equation (4.9) cannot, in general, be written in a closed form, some 
kind of approximation has to be made to determine the operator 0{t, s, z). The Markov 
approximation emerges for a vanishing environment correlation time, a{t, s) = 6{t — s). In 
this case, from (3.4) 0(t,z) = |L, and Eq. (4.9) reduces to the Markov Lindblad master 
equation, 

j^pt = -t[H,p,] + Lp,L^ - \{LL\p,} (4.10) 

Another interesting case is when the dependence of the operator 0{t,z) on the noise Zt 
is negligible, that is, 0(t,z) ^ Oo(t). Recall from (3.16) that this is indeed the case 
when the dynamics is not far from Markov or the driving noise is very small. Under this 
approximation, the master equation takes the following simple form: 

j^pt = -t[H,pt] + [L,ptOoit)^] + [Oo(t)p,,Lt] (4.11) 

The notation Oo{t) is same as before(see (3. 16), (2. 13)). The master equation (4.11) will 
serve as a good approximation to the exact non-Markovian master equation (4.9) in many 
situations of interests. In particular, if the operator 0(t, s, z) is independent of noise Z(, 
then Eq. (4.11) becomes exact. Interestingly, there are many physically relevant examples 
that satisfy this condition [1]. 

More importantly for this paper, this condition is always satisfied in the first-order 
perturbative approximation (3.4) developed in section 3.1. Then the master equation 
(4.9) takes the following form: 

j^pt = - t[H, pt] + {go{t) + 9;{t))LptL^ - 9o{t)L^Lpt - gl{t)ptL'^L 

+ tg,{t)\L\[H,L]p,]-tgl{t)\p,[L\H],L] 

+ g2{t)[L\[L\L\Lpt]+gl{t)[ptL'^[L\L\,L\ (4.12) 

This master equation is the main result of this subsection. It provides a systematic 
evolution for first-order non-Markovian systems. Hence it could be called the "post- 
Markov" master equation. As for the first-order QSD equation (3.8), the second and 
third line are smaller by a factor oot or Tt compared to the first line (recall that u is the 
typical 'system' frequency determined by //, F is a typical 'system' relaxation time scale 
determined by L^L, and r is the environmental correlation time). 

Note that this "post-Markov" equation in general remains non-Markovian even when 
giit) = 0,g2(t) = 0, because of the go(t) term. However, for long time go(t) tends to a 
constant. 

Equations (3.8) and (4.12) will be applied to some examples in Sections V. In addition. 
Section VI presents a perturbation analysis of the quantum Brownian motion model. 

Finally, it shoud be noted that we have not touched issues such as mathematical con- 
ditions for the convergence of the expansions (3.1) and (3.16). Also, we are not able to 
prove, in full generality, that Eq. (4.12) always yields a positive density operator. How- 
ever, in this paper, we shall illustrate in several examples that the resulting approximate 
QSD and master equations around the Markov limit are mathematically consistent. We 
will come back to these issues in future discussions. 
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5 Examples and Applications 

The perturbative approach developed in the previous sections allows to apply first order 
non-Markovian QSD to any open quantum system once the Hamiltonian of the system if, 
the Lindblad operator L and the environment correlation function a{t^s) are specified. 
All of these are determined by the physical model itself, as illustrated in this section using 
some typical models. For simplicity, we assume that the complex process zt entering the 
non-Markovian QSD equation (3.8) has a Lorentzian spectrum, i.e. is of the Ornstein- 
Uhlenbeck type with the correlation function a(t,s) = |e~'''l*~*l, where 7"^ = r is the 
environmental correlation time, unless otherwise stated. 

5.1 Dissipative model 

In this subsection, we consider a dissipative two-level model characterized by 

i/=|(7„ L = A(7_. (5.1) 

Since this model can be solved exactly [1], we are able to compare the perturbation ap- 
proach with the exact non-Markovian QSD and master equations. Note that the model 
defines two 'system' time scales through the parameters by 00 (oscillation) and A^ (damp- 
ing). Here we assume that they are of the same order of magnitude. 

The first-order non-Markovian QSD equation can be obtained from (3.8): 

—Jjt = - i-a:,i)t + A(o-_ - {a_)t)iptZt 
at 2 

- (\^g^{t) + i\^uog^{t) + AV(t)) {a+a_ 

- {a+)t(T^ - {a+a^)t + {a+)t{a^)t) ^t (5.2) 

and the first-order non-Markovian master equation can be obtained from (4.12): 

■^Pt = - i-[(Tz,pt] + X^go{t) {2a_pta+ - {a+a_, pt}) 

- i\'^ujgi{t)[a+a_, pt] - AV(^){o-+o"-,Pt} + 2AV(t)o-_piO-+ (5.3) 

where gi{t) gives the time-dependent frequency shift. Thus the master equation is of 
the Lindblad form with time-dependent coefficients. As seen in the next subsection, this 
property can not be regarded as a generic feature of a non-Markovian master equation. 
Note that the first-order master equation (5.3) respects the hermiticity, normalization 
and positivity for any initial states and time scales. We can easily identify the first order 
non-Markovian corrections terms of the order ui/'-f and A^/7 in equations (5.2) and (5.3). 
We expect these equations to be a good approximation for the exact solution as long as 
terms of the order (0^/7)^ and (A^/7)^ are negligible. 

In Fig. 1, the average of < <? > for uj = X = 1 and 7 = 10 are plotted. The results 
given by perturbation QSD equation over 2000 realizations (solid curve) is in remarkable 
agreement with the exact master equation (dotted curve). 
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To illustrate the limits of the Markov approximation Fig. 2 presents the ensemble 
average (ax) for the first-order QSD (solid curve) for the same parameters as Fig. 1 except 
for the memory time 7 = 1, and compares this with the Markov master equation (dotted 
curve) and the exact master equation (dashed curve). Clearly, the ensemble average of 
(Jx over 1000 trajectories still gives a good approximation to the exact master equation. 
The result is fully in accordance with our expectation as for relatively long memory times 
the Markov approximation is no longer valid. It should be noted that, in general, the 
accuracy of the first-order QSD is also limited to relatively short memory times, but 
not as severly as the Markov approximation. Then the higher order approximations or 
alternative expansion such as (3.16) should be used. 

5.2 Two-level model 

Let us consider a driven two-level atomic system interacting with a dissipative environ- 
ment. The Hamiltonian of the system, H , and the Lindblad operator, L, which represents 
the infiuence of the environment are given by 

H = |a^, L = \a, (5.4) 

respectively, where the parameter A is a coupling constant. For this model, it can be 
shown that the expansion (3.16) will not terminate at any finite order. The application of 
the perturbative approach is thus necessary. The first-order non-Markovian QSD equation 
can be readily obtained from (3.8): 

at 2 

+ A^5fo(^) [Wz)t(^z - Wz)^) i^t 

- uX'^giit) {icTx + {crz)t<yy - i{<yx)t - {(^z)My)t) i't (5.5) 

where the coefficients go(t),gi(t) are given by (3.11) and (3.12), respectively. The first 
two lines in the above equation are expected from the Markov QSD picture. The third 
line represents the non-Markovian correction and is smaller by a factor ur = ui/'j. 

Similarly, the first-order non-Markovian master equation can be obtained directly from 

(4.12): 

-nPt = - i-[(7x,pt] + 2\'^go{t)a,ptaz - 2\'^go{t)pt 

- i\^ujgi{t)[(rx,pt]- \'^ujgi{t)azpt(ry- \'^ujgi{t)ayptaz (5.6) 

There are some new features about the master equation (5.6). First, it is obviously not in 
the Lindblad form due to the presence of the cross term azpcy and its conjugate. Second, 
the master equation derived in this way naturally preserves the hermiticity, trace and 
positivity. The preservation of trace and hermiticity is obvious. It is known that positivity 
of any two dimensional density matrix is equivalent to the condition 1 1 < a > 1 1 < 1, where 
< a >= Tr(cTp) is the Bloch vector [29, 30]. In Fig. 3 we have plotted the norms of the 
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Bloch vector using the time-dependent master equation (5.6) (solid curve) and the long- 
time limit master equation(LME) (dotted curve), in which the coefficients of the master 
equation become constant: go{t) = l/2,gi(t) = I/27. Clearly, LME loses positivity for 
some initial states at short time scales, whereas, the time-dependent master equation (5.6) 
preserves postivity at all times. Note that this simple model is the 2-level analog of the 
Caldeira-Leggett master equation studied in Section VI below. 

We also solved numerically the ffist-order QSD equation (5.5). The average of a 
obtained through many realizations of (5.5) (solid line) and through the ffist-order master 
equation (5.6) (dotted curve) are plotted in Fig. 4. Taking the ensemble mean over 500 
realizations we see from Fig. 4 that the ffist order QSD equation is in good agreement 
with the ffist-order master equation, for the short memory time(7 = 10). 

6 Quantum Brownian Motion: Perturbative analysis 

The transition from non-Markovian to Markov processes is an outstanding problem. It 
is debated how to take the correct Markov limit for a non-Markovian process. Certain 
approximations of the exact dynamics can lead to master equations with bad properties 
such as non-positivity. A notorious example is the Caldeira-Leggett master equation 
[26, 31, 32] which may violate positivity of the density operator at short time scale [33, 
34, 35, 36]. Consequently, it is impossible to simulate friction a la Caldeira-Leggett 
with stochastic Schrodinger equations. The aim of this section is to apply the time- 
dependent perturbation approach for master equation developed in the previous sections 
to the Quantum Brownian Motion (QBM) model [26, 31, 32, 37]. In particular, we shall 
show that our ffist-order non-Markovian master equation recovers the Caldeira-Leggett 
master equation in the Fokker-Planck and long-time (t ^ r) limit. The Hamiltonian of 
the system and the Lindblad operator are as follows: 

H = ^p' + V{q), L = q (6.1) 

where we choose a unit mass particle moving in a general potential V{q). For the sake of 
simplicity, we consider the case of the Ohmic heat bath, I {00) ~ uj. The bath correlation 
function is then given by 

Tj [^ / UJ \ 

a{t, s) = — dcuLU coth(— — ) cos[cj(t — s)] — i sm[uj{t — s)] (6.2) 

TT JO V 2kT } 

where A is the cut-off frequency of the bath which characterizes the correlation time 
r = A^^ and r] is the friction coefficient. 
From (3.4), we get 

0{t) = g^{t)q-gx{t)v (6.3) 

where the coefficients gi^{t\g\{t) are defined as before (3. 5), (3. 6). 

The zeroth-order master equation can be obtained from (4.12) by setting g\{i) = 
g^it) = 0: 

—pt = -i[H,pt] - goR{t)[q, [q,pt]] - igoi{t)[q^, pt]- (6.4) 
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This master equation preserves positivity for all times, regardless of the initial states, as 
it is of the standard Lindblad form with time dependent coefficients. (This is of course 
not a generic feature for non-Markovian master equations). However, Eq. (6.4) does 
not take the energy dissipation into account. More relevant is therefore the first order 
approximation. The master equation in this case can be obtained from (4.12) 

—pt = - i[H,pt] - gm{t)[q, [q, pt]] - Woi{t)[q^, pt] 

+ 9iR{t) [q, [P, pt]] + igiiit) [q, {p, pt}] (6.5) 

where the coefficients giR{t), Quit) {i = 0, 1) are the real and imaginary parts of gi(t), {i = 
0, 1), respectively. The coefficient gonit) induces diffusion and the decoherence in positon 
q while goi(t) gives rise to a time- dependent frequency shift. The coefficient ginit) is 
responsible for further diffusion, and the last coefficient gii{t) gives the friction. All of 
these time-dependent coefficients vanish at t = due to the assumption that initially the 
state of bath+system is factorable. In the special case when V{q) is a quadratic potential, 
it is reassuring that our non-Markovian master equations (6.4) and (6.5) coincide with 
the zeroth and, respectively, first order expansions of the exact Hu-Paz-Zhang master 
equation [32]. 

In the Ohmic case (6.2), there exists a special high temperature limit (Fokker-Planck 
limit) which results in a Markov master equation. We take the high temperature limit in 
such a way that kT ^ A. For times t ^ t = A~^ the time dependent coefficients in (6.5) 
approach constant values and we get the Markov 'Caldeira-Leggett' master equation for 
Brownian motion: 

—Pt = -i[H',pt] -i-[q,{p,pt}] -r]kT[q,[q,pt]] (6.6) 

where H' is the cutoff-dependent renormalized Hamiltonian. This is a Markov master 
equation with constant coefficients. It does not belong to the Lindblad class and it may 
violate the positivity of the density operator. We mention casually that a next order 
high temperature expansion improves this situation and replaces the Caldeira-Leggett 
equation (6.6) by a proper Markov Lindblad equation [36]. This 1/T asymptotic expansion 
has nothing to do with the perturbative approach in our present work. Note that the 
Fokker-Planck limit of the zeroth-order non-Markovian equation (6.4) does not contain 
the dissipative (friction) term on the r.h.s. so it is a Lindblad master equation. 

It is instructive to look at the non-positivity of the Caldeira-Leggett master equation 
(6.6) from the QSD point of view. It is clear from the derivation that the QSD master 
equation (6.5) differs from the standard Caldeira-Legget master equation (6.6) for short 
times of the order of the environmental correlation time. During this short time, an 
arbitrary initial condition, which might lead to positivity violation when propagated with 
the non-Lindblad master equation (6.6), evolves towards an effective, modified 'initial' 
density operator for the long time master equation (6.6) [35]. 

Our QSD master equation (6.5) is also a non-Lindblad equation but with time- 
dependent coefficients. As in the case of the spin model in section 5.2, their time de- 
pendence can assure the preservation of the state's positivity. In the master equation 
(6.5), the coefficient gii{t) of the dissipative term is zero at t = and its time derivative 

15 



vanishes, too. The diffusion coefficient goR(t) also vanishes but its initial derivative is 
positive. Thus the initial phase of the evolution is dominated by diffusion. This mech- 
anism may, as is well Icnown in the exact model of Ref.[32], guarantee the positivity of 
the density matrix at short times as well as at later times when the dissipation enters. In 
contrast, in the Caldeira-Leggett master equation (6.6) the constant dissipative term will 
immediately violate the positivity of a distinguished class of initial density matrices. 

In summary, we have presented the zeroth-order master equation (6.4) and the first- 
order non-Markovian master equations (6.5) based on QBM model. After an initial 
'slip' time, of the order of the environmental correlation time, we recover the standard 
QBM master equation. We note that both decoherent histories and environment-induced 
decoherence are discussed using the QBM model, but mainly in the Markov regimes 
[15, 38, 39, 40]. It would be interesting to study these approaches with non-Markovian 
master equation like (6.5). We shall discuss these topics elsewhere. 

7 Conclusions 

Non-Markovian QSD offers a brand new avenue to explore non-Markovian dynamics of 
open quantum systems. Such situations appear in a variety of practical problems, like e.g. 
materials with photonic bandgaps or output coupling from a Bose-Einstein condensate. 
In this article we present a systematic perturbation approach to non-Markovian QSD. 

Our perturbation approach makes non-Markovian QSD more amenable to computer 
simulations. In particular, a detailed analysis of first-order "post-Markov" QSD equations 
and the corresponding "post-Markov" master equations are presented in Sections III and 
IV. It is noteworth that these equations depend only on the system Hamiltonian, the 
Lindblad operator and the environment correlation function. The equations can thus be 
read off directly from the total system+environment Hamiltonian. We have illustrated 
the perturbation approach with some typical examples. 

In the Markov regime, it is well-known that each Lindblad master equation can be 
unraveled by either continuous or jump trajectories which decompose the density matrix 
into pure states at all times. The reverse is also true, each stochastic unraveling uniquely 
yields a positive density matrix. In the present paper, we have shown that this corre- 
spondence is even more fruitful in the non-Markovian regime. We show explicitly how 
the non-Markovian QSD equation gives rise to the corresponding non-Markovian master 
equation. As the most important application, we have shown that each perturbative QSD 
equation naturally gives rise to a perturbative master equation. We have shown numer- 
ically that the resulting master equation naturally respects the properties of hermiticity, 
normalization and, more importantly, positivity. 

Admittedly, many issues remain to be solved in this subject. In this paper we have 
exclusively discussed the first-order "post-Markov" perturbation theory for QSD without 
touching the perturbative QSD based on the functional expansion (3.16). It is important 
to note that these two expansions (3.1) and (3.16) are of rather different physical mean- 
ing. The former expansion, on which we concentrate in this paper is an expansion in the 
environmental correlation time, whereas the latter is the expansion for the 'small noise'. 
Clearly, the comparison of these two expansions will be interesting. Another important 
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project in the next step is to apply the non-Markovian QSD to some reahstic physical 
problems such as non-Markovian atom-field interaction and in particular, the superradi- 
ance near a photonic band gap, in which the non-Markovian interaction is essential (e.g., 
see [5]). Also, it is known in the Markov regime [41, 42] that localization of quantum tra- 
jectories - typically in phase space - is of great significance in accelerating the numerical 
simulations. Therefore, investigations into localization in non-Markovian QSD would be 
useful in both theoretical and practical respects. 
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A Perturbation expansion of the operator 0(t, s, z 

Let us consider the following expansion of the operator 0{t, s, z) 

d{t,s,z) = do{t,s) 

ft ^ 
Oi{t, s,v)zydv 



10 
+ d2(t,S,Vi,V2)z^^Z^^dVidV2 

Jo Jo 

+ .... 

+ ... dn{t, S,Vi, ...,Vn)Zv^...Zv„dVi...dVn 

Jo Jo 
+ ... (A.l) 

where the operators 0„(t, s,vi, ...,Vn) are independent of the noise z and are symmetric 
in their n last variables (eg 02(t,s,vi,V2) = 02(t,s,V2,vi)). The initial condition is 
d{t,t,z) = L. 

Accordingly, we get 

- dtOo{t, .sjtpt 



rt 



dt 5z<, 

+ (di{t,s,t)zt+ dtdi{t,s,v)zydvjtlJt 

+ (2 d{t,S,t,V2)ZtZy2dV2+ / dtd2{t,S,Vi,V2)Zy^Z^2dVidV2J1pt 



+ .... 
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+ (n ... j d{t,S,t,V2,...,Vn)ZtZ,a2...Zy^dV2...dVn 

dtdn{t,S,Vi,...,Vn)z^^...Zy„dVi...dVn] i^t + ■■■ 



t rt 

Jo 



and 



+ d{t,s,z)j^iJt (A.2) 



^T+'^t = i-tH + Lzt)—- 
oz. dt oz. 



- ^o,.)^ 



- L^ (Oi{t, s) + 2 I 02{t, s, V2)z^^dv2 + ... 

+ n ... On{t,S,V2,...,Vn)z^2---Zv„dV2...dVnj'ipt + ■■■ (A.3) 



where 0„(t, vi, ..., Vn) = Jo a{t, s)On{t, s, vi, ..., Vn)ds. 
Consequently, from the consistency condition 



d 6 6 d 

dt 6zs 6zs dt 

one obtains the following hierarchy of equations: 

d - 

— On{t,S,Vi,...,Vn) = [-iH,On{t,S,Vi,...,Vn)] 

1 " - 

- — J2 II[^"^C'fc(t,wp„(i),...,wp„(fc)),0„_fe(t,s,t;p„(fc+i),...,'yp„(„))] 

^' Pn£S„k=0 

- {n + l)L^On+lit,S,Vi,...,Vn) (A.5) 

with initial conditions: 

doit,t) = L (A.6) 

On(t,t,vi, ...,Vn) = foralln>l (A. 7) 

nOnit, S, t, V2, ..., Vn) = [L, On-lit, S, V2, ..., Vn)] (A.8) 

Where S'„ is the permutation group and P„ is the permuation operators acting on the 
indices Vi,V2, ...,Vn. 

Of particular interest is n = 0, we get 

^Oo(t,s) = [-zH,do{t,s)] 



- [LtOo(t),Oo(t,s)] 

- L^Oi{t,s) (A.9) 



From this the derivative of Oq can be easily worked out. 
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From simplicity, we assume here exponentially decaying correlations: 

a(t,s) = ^e-Tl*-^l (A.IO) 

The evolution equations for the 0„(t, fi, ...,f„) read: 
d - 7 - 

— On{t,Vi,...,Vn) = -On{t,t,Vi,...,Vn)~jOn{t,Vi,...,Vn) 

+ [-iH,On(t,Vi,...,Vn)] 

1 " - 

- — Yl J2i^^0kit,Vp^(i),...,Vp^(^k)),0n-kit,Vp^(^k+l),-;Vp^(^n))] 

^- PnGSn k=0 



- (n + l)LM a{t,s)On+i{t,s,Vi,...,Vn)ds (A.U) 

Jo 

where nO„(t, t, V2, ..., Vn) = [L, On_i(t, V2, ..., fn)] for n > 1 and O„(0, t;i, ..., t;„) = for all 
n. 

The Eqs. (A. 5) and (A. 11) are very useful in the determination of the operator 

dit,s,z). 

B Second order QSD equation 

In this appendix, we present the second order non-Markovian QSD equation. 

By using the functional expansion of 0{t,s,z) and the consistency condition (See 
Appendix A), we can work out expansion of the operator 0{t,s,z) at point s to any 
desired order. In what follows, for simplicity, we only give the second order expansion of 
the operator Oo{t, s), which contains no nonlocal noise z. 

do{s,s) = L (B.l) 

9i6o(s,s) = -t[H,L]-g,is)[L^,L]L (B.2) 

dfdois,s) = -[H,[H,L]]+tgois)[H,[L\L]L]-ais,s)[L\L]L 
+tgo{s) [[L^[H, L],L]+ gl{s) [L^L^ L]L, L] 
+tgois)[L^L, [H,L]]+g^,is)[L^L, [LIL]L] (B.3) 

Note that all derivatives above are kinds of approximations, in particular, the second 
order derivative might contain more terms. Taking the first three terms of the expansion 
(3.1), and making the approximation 0{t, s, z) ~ Oo(t, s), one obtains 

Oo(t) = g^{t)L-ig^{t)[H,L\-g2{t)[L\L\L 

-g:,{t)[H,[H,L\]-g^(t)[L\L\L 

+ig,{t) {[H, [L\ L]L] + [L^[H, L],L] + [L^L, [H, L]]) 

+ge{t) ([Lt[Lt,L]L,L] + [L^L, [L\L]L]) (B.4) 



where the coefficients are as follows: 



go{t) = f a{t,s)ds (B.5) 

Jo 
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92{t) 

9A{t) 
95{t) 
9^{t) 



a{t, s){t — s)ds 



t rs 



a(t, s)a{s — u)(t — s)dsdu 



JO 

1 /■* 

- / Q;(t, s){t — s^ds 

2 JO 

1 /■* 

2 JO 

- / / a{t,s)a{s — u){t — s^duds 
2 Jo JQ 

/ / a{t,s)a{s,u)a{s,v){t — sfdvduds 

JO JO 



(B.6) 
(B.7) 
(B.8) 
(B.9) 
(B.IO) 
(B.ll) 



Then the second order QSD equations can be obtained by substituting (B.4) into (2.12). 
Notice that, in principle, we could obtain any order approximate QSD equations by di- 
rectly using the consistency condiiton and the functional expansion of 0(t, s, z). 

C Derivation of the relations (4.3) and (4.4) 

In this appendix, we shall prove the relations (4.3) and (4.4). We take (4.3) for instance: 



M[PtZt] = / dsM[ztz:]M 



6z* 



(C.l) 



Suppose the complex Gaussian measure takes the form: 



P{z)dfi = N exp 



dfi 



(C.2) 



da / dTZ*ZrP{cr, r) 

here N is the normalization constant, and /3(o", r) is a kernel reciprocal to the correlation 
function a(A,r), which is defined by 

M[z;zs] = a{t, s) (C.3) 

Note that the correlation function a(t, s) satisfies a{t, s) = a{s,t)*. We then have follow- 
ing relation: 

f a{t,Typ{T,s)dr = 5{t-s) 
Now, we consider the right hand side of (C.l): 



(C.4) 



fdsM[ztZ*]M 



5P^ 

Szt 



/f SP 

dfi dsa{t,sy-^P{z) 

-N f di2 f dsa{t,syPt—P{z) 



(C.5) 



Here, integration by parts has been used from the first line to the second line. Note that 
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6z* 
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(C.6) 



Inserting (C.6) into (C.5), changing the integration order J ds and J dr, and using the 
relation (C.4), we obtain (C.l). The relation (4.3) can be obtained by taking the hermitian 
conjugate of (C.l). 
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Captions of Figures 

Fig. 1. Ensemble average of the Bloch vector < a > over 2000 trajectories of first- 
order QSD (solid curve), and by exact master equation (dotted curve) with Hamiltonian 

H = {ijo/2)az., L = A(T_ and a{t,s) = |e~'''l*~'*l. Here we choose u = \ = 1 and 7 = 10. 
The initial state is chosen as \iPq) = |— ) + i\+). 



Fig. 2. Ensmble average oi < cTx > over 1000 realizations by using the first-order non- 
Markovian QSD (solid curve) for the same model as Fig. 3. Here, '-/ = l,uj = X = 1 and 
iV'o) = l~) + 3|+). The dashed curve is the exact mster equation for the same choice of 
parameters, and the dotted curve is the master equation in Markov limit . 



Fig. 3. Illustrtion of the norm of Bloch vector of two-level system with H = {uj/2)crx, L = 
Ao"z and expontially decaying correlation function a(t,s) = ^e"'^'*"*'. The initial state is 
chosen as the excited state l^^o) = !+)• The parameter are chosen as u; = A = 1,7 = 1/2. 
The solid curve represnts the norm of the Bloch vector by master equation. And the dotted 
line represnts the norm of Bloch vector from long-time limit master equation (LME). We 
see II < (T > II > 1 for the Bloch vector by LME for the chosen initial state at short times. 
Accordingly, LME loses the positivity at short times for some initial states. 



Fig. 4. Ensemble average of < c? > over 500 realizations (solid curve) for the same model 
as Fig. 3. Here we choose cj = A = 1,7 = 10 and the initial state \^o) = |— ) -|- v^|-|-). 
The dotted curve is an average by the first-order master equation for the same choice of 
parameters. 
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